home *** CD-ROM | disk | FTP | other *** search
/ Amiga Format CD 46 / Amiga Format CD46 (1999-10-20)(Future Publishing)(GB)[!][issue 1999-12].iso / -in_the_mag- / reader_requests / scilab / tests / matric.dia.ref < prev    next >
Text File  |  1999-09-16  |  1KB  |  80 lines
  1.  
  2. Leps=1.d-5;
  3.  
  4. //tests riccati lyapounov sylvester
  5.  
  6. //
  7.  
  8. //   lyapunov
  9.  
  10. a=rand(3,3);x=rand(3,3);x=x*x';c=a'*x+x*a;
  11.  
  12. if norm(lyap(a,c,'cont')-x) > Leps then bugmes();quit;end
  13.  
  14. c=a'*x*a-x;
  15.  
  16. if norm(lyap(a,c,'disc')-x) > Leps then bugmes();quit;end
  17.  
  18. //
  19.  
  20. //      sylvester
  21.  
  22. x=rand(3,2);a=rand(3,3);b=rand(2,2);c=a*x+x*b;
  23.  
  24. if norm(sylv(a,b,c,'cont')-x) > Leps then bugmes();quit;end
  25.  
  26. //      riccati
  27.  
  28. a=rand(3,3);
  29.  
  30. b=rand(3,2);
  31.  
  32. c=rand(3,3);c= c*c';r=rand(2,2);r=r*r'+eye;b=b*inv(r)*b';
  33.  
  34. x=ricc(a,b,c,'cont');
  35.  
  36. if norm(a'*x+x*a-x*b*x+c ) > Leps then bugmes();quit;end
  37.  
  38. aa=[a -b;-c -a'];
  39.  
  40. [xx,d]=gschur(eye(6,6),aa,'cont');
  41.  
  42. xx=xx(:,1:d);
  43.  
  44. if norm(xx(4:6,:)/xx(1:3,:)-x) > Leps then bugmes();quit;end
  45.  
  46. [xx,d]=schur(aa,'cont');xx=xx(:,1:d);
  47.  
  48. if norm(xx(4:6,:)/xx(1:3,:)-x) > Leps then bugmes();quit;end
  49.  
  50. //       discrete time case
  51.  
  52. f=a;b=rand(3,2);g1=b;g2=r;g=g1/g2*g1';h=c;x=ricc(f,g,h,'disc');
  53.  
  54. if norm(f'*x*f-(f'*x*g1/(g2+g1'*x*g1))*(g1'*x*f)+h-x) > Leps then bugmes();quit;end
  55.  
  56. aa=[eye(3,3) g;0*ones(3,3) f'];
  57.  
  58. bb=[f 0*ones(3,3);-h eye(3,3)];
  59.  
  60. [xx,d]=gschur(bb,aa,'disc');xx=xx(:,1:d);
  61.  
  62. if norm(xx(4:6,:)/xx(1:3,:)-x) > Leps then bugmes();quit;end
  63.  
  64. fi=inv(f);
  65.  
  66. hami=[fi fi*g;h*fi f'+h*fi*g];
  67.  
  68. [xx,d]=schur(hami,'d');xx=xx(:,1:d);
  69.  
  70. fit=inv(f');
  71.  
  72. ham=[f+g*fit*h -g*fit;-fit*h fit];
  73.  
  74. [uu,d]=schur(ham,'d');
  75.  
  76. uu=uu(:,1:d);
  77.  
  78. if norm(uu(4:6,:)/uu(1:3,:)-x) > Leps then bugmes();quit;end
  79.  
  80.